

run scripts/start_RDShell.m

%% Import data
opts = spreadsheetImportOptions("NumVariables", 7);
addpath('figures')

% Specify sheet and range
opts.Sheet     = "ForMatlabExport";
opts.DataRange = "A1:G3933";

% Specify column names and types
opts.VariableNames = ["DateVec", "RDPrice", "ShellPrice", "RDVolume", "ShellVolume", "RDSharesOut", "ShellSharesOut"];
opts.VariableTypes = ["datetime", "double", "double", "double", "double", "double", "double"];

% Specify variable properties
opts = setvaropts(opts, "DateVec", "InputFormat", "");

% Import
tbl = readtable("raw_data\RD Shell data.xls", opts, "UseExcel", false);

%% Convert to output type
DateVec        = tbl.DateVec;
RDPrice        = tbl.RDPrice;
ShellPrice     = tbl.ShellPrice;
RDVolume       = tbl.RDVolume;
ShellVolume    = tbl.ShellVolume;
RDSharesOut    = tbl.RDSharesOut;
ShellSharesOut = tbl.ShellSharesOut;

clear opts tbl
FinalShareCounts = [2126646,9730124];

DateVec        = DateVec(1:3933,:);
RDPrice        = RDPrice(1:3933,:);
RDSharesOut    = RDSharesOut(1:3933,:);
RDVolume       = RDVolume(1:3933,:);
ShellPrice     = ShellPrice(1:3933,:);
ShellSharesOut = ShellSharesOut(1:3933,:);
ShellVolume    = ShellVolume(1:3933,:);

meanDVs = zeros(3933,2);
counter = 0;
while (counter<3933)
   counter = counter + 1; 
   if(counter<261)
       fracDVs(counter,:) = [0 0];
   else
       %Multiply by 10 because of scaling: prices are rescaled up by 100
       %and volumes are scaled down by 1000 so 1000/100 = 10 for true value
       fracDVs(counter,:) = [10*mean(RDVolume(counter-260:counter,1).*(RDPrice(counter-260:counter,1)).*(FinalShareCounts(1,1)./RDSharesOut(counter-260:counter,1))) 10*mean(ShellVolume(counter-260:counter,1).*(ShellPrice(counter-260:counter,1)).*(FinalShareCounts(1,2)./ShellSharesOut(counter-260:counter,1)))];
   end
end

arbSizes = zeros(3933-260,1);
totalWelfare = zeros(3933-260,1);
fnEvals = zeros(3933-260,1);

counter = 260;
while(counter<3933)
   counter = counter + 1;
   Values = @(m) (RDPrice(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))-(6.86300681918852*ShellPrice(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m))));
   Values2 = @(m) (RDPrice(counter,1).*(1+(100*(1000000*m/fracDVs(counter,1))*0.0000) + (((100*abs(1000000*m)./fracDVs(counter,1)).^0.5)*0.000889.*sign(m))))./((6.86300681918852*ShellPrice(counter,1)*(1-(100*(1000000*m/fracDVs(counter,2))*0.0000) - (((100*abs(1000000*m)/fracDVs(counter,2)).^0.5)*0.000889.*sign(m)))))-1;
   
   options = optimoptions('fsolve','FunctionTolerance',1e-12,'OptimalityTolerance',1e-12,'StepTolerance',1e-6,'MaxFunctionEvaluations',10000);
   [m1,fnVal] = fsolve(Values,1,options);
   grid = sign(m1)*(0:10:abs(m1));
   evals2 = zeros(size(grid,2),1);
   counter2 = 0;
   while(counter2<size(grid,2))
       counter2 = counter2+1;
       evals2(counter2,1) = Values2(grid(1,counter2));
   end
   arbSizes(counter-260,1) = m1;
   %evals = Values(grid);
   totalWelfare(counter-260,1) = sum(evals2)*10000000;
   fnEvals(counter-260,1)=fnVal;
end

ratios = RDPrice./(6.86300681918852*ShellPrice);
ratios = ratios(261:end,:);

save('intermediate\workspace_RDShell');

%% Create figures
load('intermediate\workspace_RDShell');

fig_RDShell_welfare(DateVec, ratios, totalWelfare, optfig)
fig_RDShell_arbsize(DateVec, ratios, arbSizes,     optfig)
